General Update
The Sea Level Rise pressure layer was updated with new data for the year 2013. Previous data went through 2012. In addition, the rescaling of this pressure was recalculated to use the 99.99th percentile of raw data as the maximum. All values at or above the 99.99th percentile are given a value of 1.
Data for mean sea level rise (in mm) between January 1993 and June 2014 was downloaded on 1/12/15 from AVISO.
1. Create raster from original NetCDF file
# The following 3 lines were done once to create raster from downloaded NetCDF.
# No longer need to run these
#library(ncdf4)
#r <- raster('raw/MSL_Map_MERGED_Global_IB_RWT_NoGIA_Adjust.nc')
#writeRaster(r,'tmp/MSL_Map_MERGED_Global_IB_RWT_NoGIA_Adjust.tif')
# Read in raw data
r <- raster('tmp/MSL_Map_MERGED_Global_IB_RWT_NoGIA_Adjust.tif')
plot(r,col=cols,main='Mean Annual Sea Level Rise (mm)\n1993-2014')
2. Reproject raster to mollweide and rotate to extent -180 to 180 degrees
r = rotate(r) #rotate raster: -180 to 180
#define initial projection. Helps avoid errors when reprojecting to mollweide
projection(r) <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
mollCRS <- CRS('+proj=moll')
slr_moll <- projectRaster(r, crs=mollCRS, over=T)
plot(slr_moll,main='Sea Level Rise - Mollweide',col=cols)
3. Multiply mean annual sea level rise by number of years in dataset to get a total change in mm between 1993-2014
# (3) Multiply annual rate by all years in dataset
# data used to create r comes from January 1993 - June 2014 (21.5 years)
# Multiply values by 21+5/12 = 21.41667
slr_moll <- slr_moll*21.41667
plot(slr_moll,col=cols,main='Aggregate sea level rise (mm)\n1993-2014')
histogram(slr_moll,main ='Distribution of values')
4. All values below zero indicate a decrease in sea level rise.
We choose to disregard these as the impact of sea level rise is focused only on the increase in sea level
# (4) Clip all negative values to 0
slr_moll[slr_moll<0]<-0
plot(slr_moll,main='SLR - negative values set equal to 0',col=cols)
histogram(slr_moll,main ='Distribution of data after removing negative values')
5. Log transform the layer
slr_moll_log = calc(slr_moll,fun=function(x){log(x+1)})
plot(slr_moll_log,main='Sea level rise \nLog Transformed',col=cols)
histogram(slr_moll_log,main='Distribution of data after log transform')
6. Resample layer to 1km resolution
Original resolution of the data is 0.25 degrees. To resample the data to 1km resolution, a template ocean raster that is at the desired resolution (~1km) and contains all cells included within OHI regions, is used.
#ocean is a raster with all land clipped out - at 1km with value of 1
ocean = raster(file.path(dir_N,'model/GL-NCEAS-Halpern2008/tmp/ocean.tif'))
plot(ocean,main='Ocean Raster at 1km resolution \nUsed for resampling and clipping')
slr_1km = resample(slr_moll_log,ocean,method='ngb')
plot(slr_1km,main='SLR resampled to 1km' ,col=cols)
7. Rescale using the reference point (99.99 quantile)
All values at or above the 99.99th quantile (6.24183381230408) are set to equal 1.
#get reference point
ref = quantile(slr_1km,prob=0.9999)
#normalize by the reference point - cap all values greater than 1 to 1
r_resc <- calc(slr_1km,fun=function(x){ifelse(x>ref,1,x/ref)})
plot(r_resc,main='Rescaled Sea Level Rise Pressure',col=cols)
8. Interpolation
Interpolating the data was done using ArcGIS. All NA cells were interpolated using nearest neighbor. The python script can be found here
r_int = raster('tmp/slr_moll_log_1km_rescaled_int.tif')
plot(r_int,main='Sea level rise interpolation with nearest neighbor',col=cols)
9. Clip out ocean & Finalize layer
Using the same ocean template raster as in step 6, all land is clipped out after interpolation.
r_final = mask(r_int,ocean)
plot(r_final,main='Sea Level Rise \n Final layer',col=cols)
10. Create interpolated cells raster
The cells that were filled in using the nearest neighbor interpolation are identified in a single raster.
interp = mask(r_final,slr_1km,inverse=TRUE)
plot(interp,col=cols,main='1km cells interpolated using Nearest Neighbor')